% This function simulates using the extended perturbation method based on
% the hybrid scheme, i.e. where the standard perturbation approximation is
% used when sufficiently accurate.

function [ySim_Hybrid,xSim_Hybrid,diagOut] = simExtendedPer(model,shocks,setupEPer)

%Number of simulations
numSim = size(shocks,2);    

%% Simulation using hybrid approach
ny              = setupEPer.ny;
nx              = setupEPer.nx;
nx1             = setupEPer.nx1;
mx              = setupEPer.mx;
xSim_Hybrid     = zeros(nx,numSim);
ySim_Hybrid     = zeros(ny,numSim);
x_cu            = zeros(nx,1);
diagSim_Hybrid  = zeros(numSim,5);
Gxx             = 1/2*reshape(model.gxx,ny,nx^2);
Gxxx            = 1/6*reshape(model.gxxx,ny,nx^3);
G4x             = 1/24*reshape(model.g4x,ny,nx^4);
Gssxx           = 6/24*reshape(model.gssxx,ny,nx^2);        %scaling = nchoosek(4,2)/24
G5x             = 1/120*reshape(model.g5x,ny,nx^5);
Gssxxx          = 10/120*reshape(model.gssxxx,ny,nx^3);     %scaling = nchoosek(5,2)/120
Gsssxx          = 10/120*reshape(model.gsssxx,ny,nx^2);     %scaling = nchoosek(5,3)/120

Hxx             = 1/2*reshape(model.hxx,nx,nx^2);
Hxxx            = 1/6*reshape(model.hxxx,nx,nx^3);
H4x             = 1/24*reshape(model.h4x,nx,nx^4);
Hssxx           = 6/24*reshape(model.hssxx,nx,nx^2);        %scaling = nchoosek(4,2)/24
H5x             = 1/120*reshape(model.h5x,nx,nx^5);
Hssxxx          = 10/120*reshape(model.hssxxx,nx,nx^3);     %scaling = nchoosek(5,2)/120
Hsssxx          = 10/120*reshape(model.hsssxx,nx,nx^2);     %scaling = nchoosek(5,3)/120

residualSum     = zeros(numSim,1);
ferror          = zeros(model.ny+setupEPer.mx,numSim);

for t=1:numSim
    tic
    xSim_Hybrid(:,t) = x_cu;
    %Given x_t: compute x1_t+1 and y_t (endogenous variables)
    %Standard perturbation given certainty equivalence
    y_cu    = model.gx*x_cu;
    x_cup   = model.hx*x_cu;
    if model.orderApp > 1
        AA_cu = kron(x_cu,x_cu);
        y_cu  = y_cu  + Gxx*AA_cu;
        x_cup = x_cup + Hxx*AA_cu;
    end
    if model.orderApp > 2
        AAA_cu  = kron(AA_cu,x_cu);
        y_cu    = y_cu  + Gxxx*AAA_cu;
        x_cup   = x_cup + Hxxx*AAA_cu;
    end
    if model.orderApp > 3
        A4_cu   = kron(AAA_cu,x_cu);
        y_cu    = y_cu  + G4x*A4_cu;
        x_cup   = x_cup + H4x*A4_cu;
    end
    if model.orderApp > 4
        A5_cu   = kron(A4_cu,x_cu);
        y_cu    = y_cu  + G5x*A5_cu;
        x_cup   = x_cup + H5x*A5_cu;
    end    
    
    % computing y_t+1=y_cup
    y_cup   = model.gx*x_cup;
    if model.orderApp > 1
        AA_cup   = kron(x_cup,x_cup);
        y_cup    = y_cup  + Gxx*AA_cup;
    end
    if model.orderApp > 2
        AAA_cup  = kron(AA_cup,x_cup);
        y_cup    = y_cup  + Gxxx*AAA_cup;
    end
    if model.orderApp > 3
        A4_cup   = kron(AAA_cup,x_cup);
        y_cup    = y_cup  + G4x*A4_cup;
    end
    if model.orderApp > 4
        A5_cup   = kron(A4_cup,x_cup);
        y_cup    = y_cup  + G5x*A5_cup;
    end    
    % The residuals
    ferror(:,t)      = DSGEforesight_N([y_cu;x_cup(1:mx,1)],x_cu,y_cup,1,setupEPer);
    residualSum(t,1) = max(abs(ferror(:,t)));
    if residualSum(t,1) < setupEPer.residualMax
        %Standard perturbation
        x1_cup = x_cup(1:nx1,1);
        diagSim = [residualSum(t,1),residualSum(t,1),0,0];
    else
        %Extended path
        [Z0,y_tN,N]           = getStartingValues(x_cu,setupEPer);
        [x1_cup,y_cu,diagSim] = solutionExtendedPath(Z0,x_cu,y_tN,N,setupEPer);
    end
    timeStep               = toc;
    diagSim_Hybrid(t,:)    = [diagSim timeStep];
    %Adding derivatives with respect to sigma to correct for uncertainty
    if model.orderApp > 1
        x_cup(1:nx1,1)   = x1_cup + 1/2*model.hss(1:nx1,1);
        y_cu             = y_cu   + 1/2*model.gss;
    end
    if model.orderApp > 2
        x_cup(1:nx1,1)   = x1_cup + 3/6*model.hssx(1:nx1,:)*x_cu + 1/6*model.hsss(1:nx1,:);
        y_cu             = y_cu   + 3/6*model.gssx*x_cu          + 1/6*model.gsss;
    end
    if model.orderApp > 3
        x_cup(1:nx1,1)   = x1_cup +  Hssxx(1:nx1,:)*AA_cu + 4/24*model.hsssx(1:nx1,:)*x_cu + 1/24*model.h4s(1:nx1,:);
        y_cu             = y_cu   +  Gssxx*AA_cu          + 4/24*model.gsssx*x_cu          + 1/24*model.g4s;
    end
    if model.orderApp > 4
        x_cup(1:nx1,1)   = x1_cup +  Hssxxx(1:nx1,:)*AAA_cu + Hsssxx(1:nx1,:)*AA_cu + 5/120*model.h4sx(1:nx1,:)*x_cu + 1/120*model.h5s(1:nx1,:);
        y_cu             = y_cu   +  Gssxxx*AAA_cu          + Gsssxx*AA_cu          + 5/120*model.g4sx*x_cu          + 1/120*model.g5s;
    end
    ySim_Hybrid(:,t)     = y_cu;
    
    %Given x_t: computing x2_t+1 (Exogenous states)
    x_cup(nx1+1:nx,1) = model.hx(nx1+1:nx,1:nx)*x_cu(1:nx,1) + model.eta(nx1+1:nx,:)*shocks(:,t);
    
    %Updating to current value
    x_cu = x_cup;
end
% Adding the level
ySim_Hybrid = ySim_Hybrid + repmat(model.g0,1,numSim);
xSim_Hybrid = xSim_Hybrid + repmat(model.h0,1,numSim);

%disp(['Average time per iteration for Hybrid = ', num2str(mean(diagSim_Hybrid(:,5)))])
diagOut.maxAbsfVal   = diagSim_Hybrid(:,1);
diagOut.maxAbsfVal_1 = diagSim_Hybrid(:,2);
diagOut.iter         = diagSim_Hybrid(:,3);
diagOut.solver       = diagSim_Hybrid(:,4);
diagOut.timePerPeriod= diagSim_Hybrid(:,5);


end


%nchoosek( v , k )